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In a recent paper [F. Vega Reyes et al., Phys. Rev. Lett. 104, 028001 (2010)] we presented 
a preliminary description of a special class of steady Couette flows in dilute granular gases. In 
all flows of this class the viscous heating is exactly balanced by inelastic cooling. This yields a 
uniform heat flux and a linear relationship between the local temperature and flow velocity. The 
class (referred to as the LTu class) includes the Fourier flow of ordinary gases and the simple 
shear flow of granular gases as special cases. In the present paper we provide further support for 
this class of Couette flows by following four different routes, two of them being theoretical (Grad's 
moment method of the Boltzmann equation and exact solution of a kinetic model) and the other two 
being computational (molecular dynamics and Monte Carlo simulations of the Boltzmann equation). 
Comparison between theory and simulations shows a very good agreement for the non-Newtonian 
rheological properties, even for quite strong inelasticity, and a good agreement for the heat flux 
coefficients in the case of Grad's method, the agreement being only qualitative in the case of the 
kinetic model. 
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I. INTRODUCTION 

The development of the kinetic theory of non- uniform 
gases, extending the results by Boltzmann [l[ and 
Maxwell Q to near-equilibrium systems, started out with 
the seminal works, in the early 20th Century, by Hilbcrt 
0, Enskog 0, Chapman 0], and Burnett [al. Their re- 
sults allow for an accurate description of non-equilibrium 
states of gases (in particular, neutral gases), in the limit 
of Newtonian hydrodynamics Q (that is, small gradi- 
ents, scaled with the typical microscopic length scale, of 
the average fields). These theoretical works have been 
recently extended to the more general frame of granular 
gases where the inter-particle collisions are inelastic 0- 
|9| . The prototypical model of a granular fluid consists of 
a system of smooth inelastic hard spheres with a constant 
coefficient of restitution a. This parameter distinguishes 
ordinary gases (a = 1) from granular gases (a < 1). 

Granular matter is certainly involved, not only in many 
industrial processes but also in biological processes 
[ill ll2T | . This explains the growing interest in the study of 
granular matter. Moreover, granular flows are also chal- 
lenging from a more fundamental point of view (13 , fl3 ] . 
For instance, in the low-density regime, the Boltzmann 
equation can be generalized to granular gases. For all 
these reasons there is currently a great interest in the 
study of granular matter and a large number of research 
works have been recently published in this field (see, for 
instance, Refs. pi lij ITT1. Il5l - fl8j and references therein). 
In particular, the Navier-Stokes (NS) constitutive hydro- 
dynamic equations for granular gases have been derived 
from the Boltzmann and Enskog equations [l9l - l27j . This 
has allowed the description of important phenomena in 
granular matter, some of which were found to persist 
with the same qualitative behavior even beyond the range 
of Newtonian hydrodynamics (basic segregation mecha- 
nisms [HI, for instance). 



Unfortunately, the ranges of interest of the physics of 
granular gases fall frequently beyond Newtonian hydro- 
dynamics since the strength of the spatial gradients is 
large in most situations of practical interest (for example, 
in steady states) , due to the coupling between inelasticity 
and gradients [lj, H3|- I n these states, a hydrodynamic 
description is still valid but with constitutive equations 
more complex than the NS ones. On the other hand, 
the derivation of these non-Newtonian equations from 
the inelastic Boltzmann equation is an extremely com- 
plex mathematical task. For this reason, one is forced 
to resort to approximate schemes (such as Grad's 13- 
moment method or the use of simplified kinetic models), 
to be tested against computer simulations such as the 
direct simulation Monte Carlo (DSMC) method HU and 
event-driven molecular dynamics (MD) simulations [30j |. 
In this context, analytical solutions of the Bhatnagar- 
Gross-Krook (BGK) model kinetic equation, and its ex- 
tension to inelastic collisions, have been found for steady 
non-linear shear flows, both for clastic [3l| and granular 
gases [32Tj35| . Comparison with numerical solutions of 
the Boltzmann equation by means of the DSMC method 
shows that this kinetic model is able to describe the gen- 
eral properties of non-linear shear flows in elastic and 
granular gases. 

One of the well-known examples of steady states is 
the simple or uniform shear flow (USF) problem [l4|, Hll ■ 
This state is characterized by a linear velocity field 
(that is, du x /dy = const), constant density n, and con- 
stant temperature T . The presence of shearing induces 
anisotropies in the pressure tensor , namely, nonzero 
shear stress P xy and normal stress differences P xx — P yy 
and Pyy — P zz . On the other hand, the heat flux van- 
ishes due to the absence of density and thermal gradi- 
ents. The steady-state condition requires that the col- 
lisional cooling (which is fixed by the mechanical prop- 
erties of the granular gas particles) is exactly balanced 
by viscous heating (which is fixed by the shearing). This 
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FIG. 1: The planar Couette flow is driven by two horizontal 
plates, separated by a distance h. Both act like sources of 
temperature and shear on a low-density granular gas filling 
the space between them. 



relationship between the shear field and dissipation sets 
the strength of the scaled velocity gradient for a given 
value of the coefficient of restitution. This implies that 
the corresponding hydrodynamic steady state is inher- 
ently non-Newtonian (that is, beyond the scope of the 
NS equations) in inelastic granular gases [28| . 

Let us consider the more complex case of a generic 
planar Couette flow problem, which is depicted in Fig. 
[TJ In this state, the temperature is in principle not uni- 
form and, consequently, a heat flux vector q coexists with 
the pressure tensor P,j [33j . In fact, the energy balance 
equation (in the steady state) reads 



dqy 
dy 



p dux_ 

XV dy ' 



(1) 



where Q is the inelastic cooling rate and d is the dimen- 
sionality of the system. The first term on the right- 
hand side is an energy sink term reflecting the dissipa- 
tion due to collisions, while the second term [note that 
sgn(P xy ) = — sgn(du x / dy)} is an energy source term due 
to viscous heating. The competition between these two 
terms determines the sign of the divergence of the heat 
flux [36| . As for the conservation equation for momen- 
tum, it implies 



P xy = const, 



Pyy = COnSt. 



(2) 



(3) 



In general, Eq. (JXJ) applies to any state that (i) is sta- 
tionary, (ii) has gradients only along the y direction, 
and (iii) has a flow velocity vector along the x direction. 
Thus, Eq. ([1]) is also valid for the familiar Fourier flow of 
ordinary gases (a = 1) as well as for the (steady-state) 
USF of granular gases (a < I). In the first case C = 
and du x /dy = 0, so the non-zero heat flux vector is uni- 
form. In the second case, there is no heat flux and, as 
said before, the condition 



c 



2 P 8v 
d nT dy 



(4) 



establishes the relationship between the inelastic cool- 
ing and the shear field. These two clearly distinct states 



share the common features of uniform heat flux and a 
local balance between inelastic cooling and viscous heat- 
ing. The interesting question is, docs there exist a whole 
class of Couette flows also sharing the same features? 
This class would include the Fourier flow of elastic gases 
and the USF of inelastic gases as special limit situations. 

The aim of this paper is to provide numerical and ana- 
lytical evidence on the existence of such a class of Couette 
flows. On the numerical side, we have solved the inelas- 
tic Boltzmann equation by means of the DSMC method 
[29j and have carried out MD simulations of dilute gran- 
ular gases. On the analytical side, we have solved this 
special class of Couette flows from a simplified model ki- 
netic equation as well as by the application of Grad's 13- 
moment method to the Boltzmann equation. A further 
theoretical support for this class has recently been found 
from an exact solution of the Boltzmann equation for in- 
elastic Maxwell models [3?| . Apart from the condition 
q = const, this class of Couette flows is macroscopically 
characterized by a uniform pressure, 



p = nT = const, 



and 



-%T 



y u>x 



A = const, 



a{a) = const, 



(5) 



(6) 



(7) 



where v oc nT 1 ! 2 is an effective (local) collision frequency. 
As a consequence of Eqs. © and ([7]), while neither u x (y) 
nor T{y) are linear, a parametric plot of T vs u x shows a 
linear relationship. For this reason, we refer to this class 
of flows as "linear T(u x )" flows, or simply, "LTu" flows. 
The slope of the linear plot T(u x ) goes from zero in the 
inelastic USF limit (constant temperature) to infinity in 
the elastic Fourier flow (zero macroscopic velocity). As 
we will see, the transport properties in the LTu class 
arc highly non-Newtonian and can be characterized by 
a generalized shear viscosity, normal stress differences, a 
generalized thermal conductivity, and a cross coefficient 
associated with the x component of the heat flux. A 
preliminary report of the LTu has been published recently 

The paper is organized as follows. In Sec. [II] we present 
the formal description at a kinetic theory level of the LTu 
flows, derive the relation between the Reynolds number 
and the relevant parameters, and define the generalized 
transport coefficients. We find in Sec. Mil two analyti- 
cal solutions of the problem introduced: in Sec. IIII Al an 
approximate analytical solution is obtained by means of 
Grad's 13-moment method, whereas in Sec. IIII B"l we find 
an exact solution of a model kinetic equation (BGK-type 
kinetic model adapted to the granular gas j39|). In Sec. 
HVI thc simulation techniques (both DSMC and MD) used 
in this work are described. Theory and simulation results 
are shown and compared in Sec. [V] Finally, in Sec. IVII 
wc give a brief summary of results and discuss them. 
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II. BOLTZMANN DESCRIPTION OF THE LTU 
FLOW 



-V ■ P = 0, 



ran 



(15) 



A. Couette flow 

Let us consider a granular fluid modeled as a gas of 
inelastic hard spheres. A constant parameter, the coeffi- 
cient of normal restitution a, accounts for the inelasticity 
in collisions. Its values range from a = (purely inelas- 
tic collision) to a = 1 (purely clastic collision). In the 
low-density regime, the one-particle velocity distribution 
function /(r, v; t) obeys the inelastic Boltzmann equation 
1,0 



(9 t + vV)/(r,v;i) = J[v|/,/], 



(8) 



where the Boltzmann collision operator J[v|/, /] is given 
by 

J[vi|/,/] = a*- 1 J dv 2 J dBQ(g-a)(g-a) 

x[a- 2 /K)/K)-/( Vl )/(v 2 )]. (9) 

Here, a is the diameter of a sphere, Q(x) is Hcaviside's 
step function, <r is a unit vector directed along the cen- 
ters of the two colliding particles, g = V! — v 2 is the 
relative velocity, and the primes on the velocities denote 
the initial values {v^,v 2 } that lead to {vi,V2} following 
a binary collision: 



a J ) (£-g)£, 



v 2 + - (1 + a x ) (ct • g)S\ 



(10) 



At a hydrodynamic level, the relevant quantities are 
the number density n, the flow velocity u, and the gran- 
ular temperature T. They are defined as moments of the 
velocity distribution as 



J dv/(v), 



u = - / dvv/(v) 



= ^ / dv^ 2 /(v) 



(11) 



(12) 



(13) 
u(r) is 



where m is the mass of a particle and V = 
the peculiar velocity. 

The Boltzmann collision operator conserves the num- 
ber of particles and the momentum, but the kinetic en- 
ergy is not conserved. The corresponding balance equa- 
tions are obtained by multiplying both sides of Eq. ([8]) 
by 1, v, v 2 , and integrating over velocity. The result is 



D t T + — (V • q + P : Vu) = -<T. 

an 



(16) 



Here, D t = dt + u • V is the material time derivative, 



Pi, = m 



is the pressure tensor, 



= _ / dvF 2 V/(v), 



is the heat flux, and 



dnT 



dvF 2 J[v|/, /] 



(17) 



(18) 



(19) 



is the cooling rate characterizing the rate of energy dis- 
sipated due to collisions. 

In the planar Couette flow the granular gas is enclosed 
between two parallel, infinite plates (normal to the y axis) 
at y = ±h/2 in relative motion along the x direction, and 
kept, in general, at different temperatures (cf. Fig. [T]). 
The resulting flow velocity is along the x axis and, from 
symmetry, it is expected that the hydrodynamic fields 
only vary in the y direction. Consequently, the velocity 
distribution function is also assumed to depend on the 
coordinate y only. Moreover, we focus on the steady 
state, so Eq. ((8j) becomes 



VyOyf = J[v|/, /]. 



(20) 



Under the above conditions, the mass conservation equa- 
tion p4[) is identically satisfied, the momentum conser- 
vation equation (TT5"|) reduces to d y Pi y ~ [cf. Eqs. @ 
and ©], while the energy balance equation (jTB")) becomes 
Eq. ([1]). It must be noted that Eqs. (P)-© are exact 
consequences of the geometry of the problem and the 
steady-state condition. Therefore, they are valid whether 
a hydrodynamic description applies or not, even near the 
walls where boundary effects are not negligible. 

Now we assume that the separation h between the walls 
is large enough (that is, it comprises a sufficient number 
of mean free paths) as to identify a bulk region where 
a hydrodynamic description is expected to apply. Here 
the term "hydrodynamics" is employed in a wide sense 
encompassing both Newtonian and non-Newtonian be- 
havior. In the context of the Boltzmann equation, a hy- 
drodynamic description is linked to a normal solution, 
namely, a special solution where all the space and time 
dependence of the velocity distribution function takes 
place via a functional dependence on the hydrodynamic 
fields [H: 



D t n + ?7.V ■ u = 0, 



(14) 



/ = /[v|n,u ! T]. 



(21) 
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B. LTu flow 



In the general Couette flow problem, the imposed ve- 
locity and temperature gradients can be controlled inde- 
pendently of the coefficient of restitution via the bound- 
ary conditions. This problem was studied by means of a 
simple kinetic model in Ref. [33j . Here, however, as said 
in the Introduction, we focus on a special class of Cou- 
ette flows. More specifically, we assume that there exists 
a normal solution of the Boltzmann equation (|20[) with a 
uniform heat flux component q y . As a consequence, the 
shear rate du x /dy is not a free parameter but it is fixed 
by the value of the coefficient of restitution [cf. Eq. Q]. 

As indicated by Eq. (f2"Tj) . we need to specify the form 
of the hydrodynamic fields in order to characterize the 
normal solution corresponding to the class of Couette 
flows with uniform heat flux. This is a non-trivial risky 
task since the proposed spatial dependence of the fields 
must be consistent with Eqs. © and ([3]) and, moreover, 
the state is expected to lie outside the realm of the NS 
regime. 

We take two basic assumptions (which have already 
been shown to be fulfilled for generic Couette granular 
flows UK). First, the exact condition is extended to 
the remaining diagonal elements of the pressure tensor, 
so that its trace is also uniform. This first assumption 
is displayed in Eq. ©. Note that in the NS description, 
P = Pyyi so Eq. © is a straightforward consequence of 
the conservation of momentum. Here, however, we as- 
sume Eq. (J5j) even though, as will be seen below, p ^ P yy . 
The second assumption is subtler and refers to the y com- 
ponent of the heat flux. According to the concept of a 
normal solution q y = q y [n, u, T] is a functional of the 
hydrodynamic fields. We assume that such a functional 
dependence has the same form as in the NS description, 
namely, q y cx (p/nT 1 / 2 )9 1 ,r. Note, however, that the 
proportionality constant is in general different from the 
NS one. Since p has already been assumed to be uni- 
form and q y = const defines the LTu state, it follows Eq. 
(p|) with v cx nT 1 / 2 . Therefore, Eqs. © and © define 
the assumed hydrodynamic profiles. The energy balance 
equation © yields Eq. ([7]), where we take into account 
that £ cx v as well as Eqs. © and (JSJ) - The constant 
a(a) is a dimensionless parameter that plays the role of 
the Knudsen number (Kn) associated with the shearing. 
As indicated by the notation, a(a) is not a free parame- 
ter but depends on the coefficient of restitution through 
Eq. . On the other hand, the constant parameter A 
defined by Eq. ([6]) is not constrained by the value of a. 
Note that A is not a dimensionless number, the corre- 
sponding Knudsen number associated with the thermal 
gradient being e = A/y/mT. 

As said in Sec. HI from Eqs. (|6|) and one obtains 



dT 

du x 



A 

a(a) 



= const. 



(22) 



This means that if the spatial coordinate y normal to 
the moving plates is eliminated between temperature and 




'"e/a„.. 



(Equilibrium) 



FIG. 2: (Color online) Each point of this diagram represents 
a Couette flow steady state. The surface defines the LTu 
class, which contains the lines representing the Fourier flow 
for ordinary gases (that is, no shearing and no inelasticity) 
and the USF for granular gases (that is, no thermal gradient). 



flow velocity the resulting profile T(u x ) is linear, thus 
justifying the acronym LTu used here to refer to this class 
of flows. 

It is interesting to get the explicit spatial dependence 
of T and u x [36[ • From Eq. ([6]) it is easy to obtain 



T(y) = T Q 



3Avq 
~2To~ 



(y - yo) 



2/3 



(23) 



where yo is an arbitrary reference point in the bulk re- 
gion, and Tq and vo are the local values of T and v, 
respectively, at y = y . Integrating Eq. (f22j) , with the 
aid of Eq. ([231) , we simply get 



Ux[y) = ~^~ T o 



3Av 



{y - vo) 



2/3 



a(«) T 
""o -j- T , 



(24) 

where uq is the local value of u x at y = yo. The expres- 
sion of the (local) thermal Knudsen number is 



<v) 



A 



\fmTo 



iAuo 



l -i/3 



(y - yo) 



(25) 



In the particular case of elastic particles (a = 1 or, 
equivalently, £ = 0), Eq. (|4]) implies a = 0, so we recover 
the Fourier flow of ordinary gases |4£j. On the other 
hand, in the absence of thermal gradients (i^O) but in 
the presence of inelastic collisions (a < 1), Eqs. (|23|) and 
(|24[) become T = To and u x = uo + a(a)vo(y — yo), that 
is, we recover the conditions of USF. For general values 
of a and A, Eqs. ©-(JT]) define a whole class of Couette 
flows with uniform q y . This manifold of Couette states 
is sketched in Fig. [2] On the LTu surface one has d y q y = 
0, while the points above (below) the surface represent 
Couette-flow states where the dominant term in Eq. ([T]) 
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is the viscous heating (inelastic cooling) one and thus 
d y q y > [d y q y < 0). For an analysis of the curvature of 
the temperature profiles within the NS domain, see Ref. 

U. 



C. Reynolds number for LTu flows 

So far, we have not needed to specify the explicit form 
of the effective collision frequency v, except for the scal- 
ing relation v cx nT 1 / 2 . Henceforth, we will adopt for v 
the conventional choice of effective collision frequency in 
shear flow problems involving ordinary gases, namely, 



(26) 



where 77^ is the NS shear viscosity of a gas of elastic 
hard spheres. With this choice, one has (in the leading 
Sonine approximation) [|| 



87r (d-l)/2 



-n\l —a 

m 



.d-l 



(27) 

(d + 2)T( d 

It is instructive to express the Reynolds number of the 
LTu flow in terms of the reduced shear rate a(a), the wall 
temperatures T±, the slab width h, and a nominal mean 
free path I. The Reynolds number Re is defined as [4l[ 

Re= mn(U + _-U-)h ^ (28) 

where n and fj^ s are characteristic values for density and 
shear viscosity, respectively. Here we take n as the aver- 
age number density and 77^ = p/v, where v is given by 
Eq. ([27]) by setting n = n and T = T_ . 

Neglecting velocity slips and temperature jumps near 
the walls, and choosing yo = —h/2 in Eqs. (|23|) and ([24]) . 
one obtains 

U + -U. = ^(T + -T_) 
3 AT 



-a(a)u(-h/2)h, (29) 
2(1 + AT) 3 / 2 - 1 ' ' v ' 

where AT = T + /T_ — 1 and, without loss of generality, 



we have assumed T + > T_ . Insertion of Eq. (|29)) into Eq. 
([28]) yields 

3 AT , . fh\ 2 , . 

Re= 2 (l + AT)3/^i a ^( v ?J ' ( 3 °) 



where i = y/T-/m/ v is the nominal mean free path. 
Upon derivation of Eq. (|30| use has been made of the re- 
lation u(-h/2)/v = n(-h/2)/n = p/nT_. Equation ([SO]) 
expresses the Reynolds number in terms of the relative 
temperature difference AT, the shear-rate Knudscn num- 
ber a(a), and the system-size Knudscn number £/h. We 
observe that Re is essentially the ratio between the shear- 
rate Knudsen number and the square of the system-size 
Knudsen number. The pre-factor depends on AT and 
ranges from 1 in the limit AT — > to in the opposite 
limit AT — > 00. 



D. Non-Newtonian transport coefficients 

As said above, the LTu flow is in general non- 
Newtonian. This can be characterized by the introduc- 
tion of generalized transport coefficients measuring the 
relationship between momentum and heat fluxes with 
the hydrodynamic gradients. First, we define a non- 
Newtonian shear viscosity coefficient 77(a) by 



P, 



-77(a) 



du x 



(31) 



Since, by dimensional analysis, 77 cx p/v 1 Eq. pip is con- 
sistent with Eqs. ©, ©, and ((TJ). Equation (JHTJ can 
be seen as a generalization of the NS constitutive equa- 
tion for the shear stress in the sense that it is assumed 
that P xy is independent of the thermal gradient A. On 
the other hand, the generalized shear viscosity coefficient 
77(a) is expected to differ from the NS shear viscosity co- 
efficient ?7ns (a) of an inelastic dilute gas [2l| . The energy 
balance equation (01 establishes a relationship between 
the reduced shear rate a(a), the generalized shear vis- 
cosity 77(a), and the cooling rate C( a ) : 



a 2 (a) 



d(*(a) 
2 if (a)' 



(32) 



where Q* = C,/v and 77* = rj/(jp/i>). 

While P xx = P yy = p in the NS regime, normal stress 
differences are expected to appear. They can be mea- 
sured though the coefficients 



P.. 



XX 

P 



.(a), 



' vv 
P 



9 v (a). 



(33) 



For d > 3, one could define a coefficient 9 Z = P zz /p but 
it is related to 8 X and 6 y by the condition 9 X + 6 y + (d — 
2)9 Z = d. The quantities 6 X and 6 y represent directional 
temperatures T x — P xx /n and T y = P yy /n (relative to 
the granular temperature T) along the x and y directions, 
respectively. 

In the case of the heat flux, the assumed scaling rela- 
tion q y cx (p/v)d y T suggests the Introduction of a gener- 
alized thermal conductivity coefficient A(a) as 



dT 



(34) 



This equation has the same form as Fourier's law, except 
that the coefficient A(a) is expected to differ from the 
corresponding NS thermal conductivity coefficient of an 
inelastic dilute gas [2l[. Moreover, while q x = in the NS 
description, here we assume the existence of a nonzero 
x component of the heat flux due to a non-Newtonian 
coupling between shearing and temperature gradient. To 
characterize this non-Newtonian effect, we introduce a 
cross coefficient 6(a) as 



q x = 6(a) 



dT 
dy ' 



(35) 
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Dimensional analysis shows that A oc p/v and <f> oc p/v, 
so that Eqs. (fMf and (|3"5j) imply that q is uniform. 

It must be borne in mind that in this section we 
have assumed the existence of Couette flows with (a) 
q y = const and (b) profiles given by Eqs. ©-([7]), but 
there is no a priori guarantee that the Boltzmann equa- 
tion ([20)1 admits such states. In the next section we will 
provide support for the existence of this LTu class by solv- 
ing Eq. (|20[) through the approximate Grad 13-moment 
method and by an exact solution of a model kinetic equa- 
tion of the inelastic Boltzmann equation. Further sup- 
port will be given by computer simulations, showing a 
good agreement with some of the theoretical results. 

III. THEORETICAL APPROACHES 



m/dvVW,/-^!^-,*). (39) 

Moreover, the collisional moments associated with the 
momentum and energy transfers are approximated by 

m J dV ViV } J[f, f\^-fav (Pij ~ pSij ) - {Pij , (40) 

j J dW 2 VJ[/,/] H- -^(huq, (41) 

where 

C = ^(l-« 2 )> (42) 



A. Grad's moment method 

In order to check the consistency of the hydrodynamic 
profiles ([5])— dTJ) , as well as of the momentum and heat 
fluxes, here we will solve the Boltzmann equation by the 
classical Grad moment method \vi ■ This in turn will 
provide explicit expressions for the generalized transport 
coefficients rj, 9i, A, and <f>. 

The idea behind Grad's moment method is to expand 
the velocity distribution function / in a complete set of 
orthogonal polynomials (generalized Hermitc polynomi- 
als), the coefficients being the corresponding velocity mo- 
ments. Next, the expansion is truncated after a certain 
order k. When this truncated expansion is substituted 
into the hierarchy of moment equations up to order k 
one gets a closed set of coupled equations. In the stan- 
dard 13-moment method the retained moments are the 
hydrodynamic fields (n, u, and T) plus the irreversible 
momentum and heat fluxes (Pij — pSij and q). More 
explicitly, 



/ -> /o 1 + 




where 



f ° = n i£f) 



m_\ d/2 ~ m V 2 /2T 



(36) 



(37) 



ft = 



d-l, 



16 + lld-3(d + 8)a 
16(d- 1) 



(1 + a). 



(43) 



(44) 



It is important to remark that, upon writing Eqs. (|40)) 
and (|4ip . nonlinear terms in P^ — p8ij and q are ne- 
glected. This is the usual implementation of Grad's 
method, although the quadratic terms are sometimes re- 
tained |43l . [4^ | . Note that the expression of the cool- 
ing rate £ provided by Grad's method and given by Eq. 
(|42|) coincides with its local-equilibrium form. The di- 
mensionless parameters (3i and P2 measure the impact of 
inelasticity on the collisional transfer of momentum and 
energy, respectively. Both coefficients reduce to unity in 
the elastic limit. 

Now, let us apply Grad's method to the Boltzmann 
equation (|2"D|) . In the geometry of the Couette flow, the 
relevant moments are n, u x , T, P xy , P xx , Pyy, q x , and 
q y . Of course, the exact balance equations ([l])-([3]) ar( 3 
recovered. The remaining five equations are obtained 
by multiplying both sides of Eq. flU by V x V y , V 2 , Vy 2 , 
V V x , and V 2 V y , integrating over velocity, and applying 



(45) 



the approximations (|38j) - (|4T|) . The results are 



d+2 



d y q x 



(filV + QPsy, 



is the local equilibrium distribution. In the three- 
dimensional case, there are 13 moments involved in Eq. 
(f36|) ; hence this method is referred to as the 13-moment 
method. In the case of a general dimensionality d the 
number of moments is d(d + 5)/2 + 1. 

In order to have a closed set of equations for n, u, T, 
Pij — pSij , and q we need to make use of Eq. (|36f to get 



in f 

- J dvViVjVkf 



d + 2 



(q*5jk + '/..'>./ + QkSij) , (38) 



d + 2 



G 



>>■:/',, 1" (46) 



d+2 



dyQy = -Pi" ( p yy - p) - CP 



2 Mm 9 



d + 4 



q y d y u x = - 



d-l 



(47) 



foqx, (48) 
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d + 2 



T {d + A 

P yy -p) + j-^q x clyu x = —fevqy, 

(49) 

Wc have made no extra assumptions in the set of equa- 
tions (|45|) — (j49]) obtained within the Grad method, apart 
from the stationarity of the system and the geometry and 
symmetry properties of the planar Couette flow. Now we 
look for hydrodynamic LTu solutions, that is, solutions 
consistent with q = const and Eqs. ©-(J?]). It is easy to 
check that Eqs. (|55 ]) -(|g9"] ) . together with Eq. ((J), indeed 
allow for such a class of solutions. First, Eqs. (|4"5j) - (|4"T)) 
become a set of algebraic equations whose solution yields 
Pxy/p, P-xx/p, and Pyy/p in terms of a and a. The re- 
duced shear rate a is subsequently obtained as a function 
of a from Eq. (|32p . Once the pressure tensor is known, 
Eqs. (|4"8|) and (|4T)|) provide q x /A and q y /A as functions 
of a for arbitrary A. The results can be conveniently 
expressed in the forms of Eqs. fl3TJl, ([3"3"]). and ([55)1 

with the following explicit expressions for the generalized 
transport coefficients: 



(/?i+C*) 5 



Pi + dC 



(50) 



(51) 



(52) 



A* = (3; 



(d - l)(d + 2)[(d + A)6 y -2}+d 2 (d + 4)(C7/3 2 ) 



(d + 2)2(d-l)/3 2 2 -2^±i) a 2 



(53) 



c/)* = (d + 4)a 



d[(d + 4)0„ - 2] + (d - l)(d + 2)77*^2 
(d + 2)2(d-l)/3 2 2 -2^) a 2 



(54) 

Here, we recall that rj* = rj/(p/v) and C* = C/^- Accord- 
ing to Eq. (|32|) , the dependence of the reduced shear rate 
a(a) on the coefficient of restitution a is 



dC 



* \2 



(55) 



In Eqs. (j53"j) and ([S"4"]l we have introduced the reduced 
coefficients A* = A/A^ s and </>* = </>/A^ g , where 



A 



_ d(d + 2) p 



NS 



2(d- 1) mi/ 



(56) 



is the NS thermal conductivity in the elastic limit. As a 
simple test, note that in the limit a — > 1 (that is, f — > 0) 
one has a — > 0, — ► 1, 0, — > 1, ?7* — >• 1, A* -> 1, and 

4>* ->• o. 

From the symmetry relation 9 X + 9 y + (d — 2)6 z = d 
and from Eqs. (|ST1) and (|52"j) it follows that 2 = 



Equations (|5C^ - (|55)) extend to arbitrary dimensionality 
d our previous results for hard spheres (d = 3) [Hj] . 

The transport coefficients (|5(3 |) -(f5"4 |) describe the non- 
Newtonian properties of the granular gas in the LTu class 
of flows in the context of Grad's solution to the Boltz- 
mann equation. These coefficients clearly contrast with 
the ones obtained in the NS description, where one has 

EMU 



?7nS 



5(7 



A NS — 



2{il 



— 1) 



(^2 



2d 



iC* A 



:?(/ 



2(d- 



— C* 
1) 



(57) 



(58) 



Upon writing Eq. (|58l) we have taken into account that 
the NS constitutive equation q = — knsVT — /insVtt. 
becomes q = — AnsVT, with Ans = k ns — ('VT')mns, 
under the condition Vp = 0. In Eqs. (|57p and (|55|) . non- 
Gaussian corrections to the homogeneous cooling state 
distribution have been neglected, in consistency with the 
Grad approximation (f36|) . Apart from Eqs. (f57| and ([58]). 
the NS description predicts #j = 1 and <fi = 0. 

Figure [3] compares the non-Newtonian coefficients 
rf (a) and A* (a) with their NS counterparts ?7Ns( a ) an d 
A NS (a) f° r hard disks (d = 2) and hard spheres (d = 3). 
It is apparent that the LTu shear viscosity clearly differs 
from the NS shear viscosity. In fact, while the latter in- 
creases with increasing inelasticity, the former presents 
the opposite behavior [28]. On the other hand, both 
thermal conductivity coefficients arc rather close to each 
other, especially in the case of hard spheres. It is inter- 
esting to remark that, while the NS heat-flux transport 
coefficients kns and /xns increase with inelasticity, the ef- 
fective NS thermal conductivity Ans = k ns — ( n /T)p^s 
decreases. This shows the importance of the coefficient 
/ins (absent in the elastic case) in granular flows beyond 
the quasi-elastic limit. 



B. BGK-type kinetic model 

Now we consider the results derived for the LTu class 
from a BGK-type kinetic model of the Boltzmann equa- 
tion [39[ . In the geometry of the Couette flow, the steady 
kinetic model becomes 



df 
dy 



^ = -/%>(/ 



f) + C d 

/o)+ 2dv" 



V/, 



(59) 



where v is the effective collision frequency defined by Eq. 
(|27|) . The parameter /3(a) is a free parameter of the 
model chosen to optimize the agreement with the Boltz- 
mann results. In terms of the variable s(y) defined as 
ds = f3v(y)dy, Eq. (|5"9")) can be rewritten as [3S 



d 







1 



l-ot + VygZ-aVy— - -CV- — ]/( S ,V) 



2' ' dV / 
/o(s,V), 
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FIG. 3: (Color online) Reduced shear viscosity (n*) and ther- 
mal conductivity (A*) for inelastic hard disks (top panel) 
and hard spheres (bottom panel), as obtained from Grad's 
13-moment method (solid lines) and from the NS equations 
(dashed lines). 



where a = a//3, C = C* IP, an d the derivative d s is taken 
at constant V = v u(s). Upon writing Eq. use 



has been made of Eq. (JT)). The hydrodynamic solution 
to Eq. (jnn]) is 



^ te -(l-iOt e -r(t)V y d Be atV y a Vx 



x/o(*,e* ft V), (61) 



where 



r(t) = i (e*& - l) . (62) 

The action of the operators e~ TVvds and e atVydv * on an 
arbitrary function g(s,\) is [33| 



e-^g( S ,V) = g(s-T^-,\) 



Vu 



(63) 



e atv v d Vjs 5 ( SjV ) = ff (s,V + art^x), (64) 
respectively. The solution (|61[) adopts the normal or hy- 
drodynamic form since its spatial dependence only occurs 
through a functional dependence on the hydrodynamic 
fields n(s), u(s), and T(s) via the local equilibrium dis- 
tribution /o. 

The objective now is two- fold. First, we want to check 
that the LTu profiles ©-([7]) are consistent with the so- 
lution ([6i"|) . Next, we will evaluate the fluxes and iden- 
tify the generalized transport coefficients defined by Eqs. 
([3T]). (|55|) . and ([55]). In order to accomplish this 

two-fold objective, it is convenient to define the general 
velocity moments 



M klMM {s)= J dVV x k W^V z k *f( S ,V). 



(65) 



Insertion of Eq. (|61j) yields 



Mi 



df / dVe- (1 -' c ) t (V a: -otV r „) fcl V i f 2 V;' C3 e- T ( t ) 1/ '' a =/ ( S ,e5C*V) 
dt / dVe-( 1+ 50*(y x -ai^) fel K fc3 V; fe3 e- ri W^ as /o(s ) V), 



(66) 



where k = ki + &2 + &3 and T"x(t) = T~(t)e af* = 2 ^1 — e aC*^ /£. It is now convenient to expand the operator 
e -ri(t)v B a s ^ go (jggj) becomes 



(') = EftlEry / d*e-^*Q* [-n(t)] h (-att-*d* I dvvX 



1=0 

fci 



1 



h=0 



K K3 fo(s,V) 



1=0 



h=0 



k+ h-l-k 3 Ck 3 



-4 



k,h,ki-td s 



n{s) 



2T(s) 



(k+h)/2 



(67) 
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where 



= U- 1/2 r(^i), t = even, 
£ | 0, I = odd, 



(68) 



.4 



3,1,2 



4a 2 (12 + 42C + 37C 2 ) 
(1 + 2C) 3 (2 + 3C) 3 ' 



(74) 



and 



A k , hM = / dte-^^hnWl^-ai) 
Jo 

In particular, A) ; o,o = 1, 

1 . 2 

^2,0,0 — =i ^2,0,1 : 



1 + C 



^2,0,: 



(i + C) 5 



2S 2 



(i + C) E 



i3,l,0 



^3,1,1 



(1 + 2C)(2 + 3C) 



2a(4 + 7Q 
(1 + 2C) 2 (2 + 3C) 5 



(69) 



(70) 



(71) 



(72) 



(73) 



.4 



3,1,3 



125 3 (4 + 70(8 + 28C + 25C 2 
(l + 20 4 (2 + 3C) 4 



(75) 



Note that because of the parity properties of the coeffi- 
cients Ci , only the terms with £ = even and h + k = even 
contribute to the summations in Eq. (|67[) . Moreover, the 
moments Mk lt k 2 ,k 3 with &3 = odd vanish. 

So far, no specific spatial dependence of density and 
temperature has been assumed. Only the linear s depen- 
dence of the flow velocity has been used. Now, we assume 
that n(s)T(s) = const and d s T(s) = const, in agreement 
with Eqs. ([5]) and ((6j), respectively. These assumptions 
imply that d^[T(s)}( k+h - 2 ^ 2 = if h > (k + h - 2)/2. 
Therefore, the summation Y^hLo can ^ e replaced by 
XX=o < '° an d Eq. (|6"T|) reduces to 



M klMM {s) = n(s) 



2T(s) 



k/2 /ci 



E 



. max(0,fc-2) 
K\ \ <^l ( ~'k+h,-C-k 3 <^k 3 



E 

h=0 



h\ 



A 



k,h,ki-l- 



( k+h 
\ 2 

k-h 
2 



ill 



mT 



AT 



(76) 



It is straightforward to check that M 0) o,o( s ) = n(s) and 
A/1.0,0 = -^0,1,0 = Ajf ,o,i = 0. This proves the consis- 
tency of the assumed density and velocity profiles in the 
LTu flow. The consistency condition for the tempera- 
ture is Af 2i o,o + Mo 2,0 + (d ~ 2) A/0,0.2 = dp/m. It can 
be checked that this condition is satisfied provided that 
the reduced shear rate a is related to the coefficient of 
restitution by 



(77) 



« 2 = |C(i + 2 - 



This result is fully equivalent to Grad's prediction (f55|) . 
except that /3\ is replaced by f$. 

Once we have proven that the BGK-typc kinetic equa- 
tion (|59|) admits an exact solution characterized by the 
LTu hydrodynamic fields, we can obtain all the veloc- 
ity moments from Eq. (|76[) . The relevant elements of 
the pressure tensor are P xx = mMififi, P yy = toA'/o,2,Oj 
and P xy = mMi^o. From them one can easily identify 
the dimensionless coefficients defined by Eqs. (f5T|) and 
(|3"3"|) . The resulting expressions coincide with Grad's re- 
sults (|5"01) - (|5"2")) . again with the replacement f3\ —> f3. 



The two non-zero components of the heat flux are 
q x = (m/2) [Af 3 ,o,o + A/1,2,0 + (d- 2)Mi j0i2 ] and q y = 
(m/2)[Mi j2 ,o + Mo,3,o + (d-2)Mo,i, 2 ]. As expected, 
they are proportional to the temperature gradient and 
this allows one to identify the generalized thermal con- 
ductivities defined in Eqs. (|3"4")l and (|3"5"|) . After some 
algebra, one gets 



A* = 



2/73 



(1 + 20(2 + 30 



2a 



1 



6a 2 12 + 42C + 37C 2 
rf+2(l + 2C) 2 (2 + 3C) 2 



rf + 2(l + 20 2 (2 + 30 2 

+ 28C + 25C 2 



18a 2 



(l + 20 2 (2 + 30 s 



(78) 



, (79) 



where we have taken into account that in the BGK model 
the NS thermal conductivity in the elastic case is not 
given by Eq. (|56|) but by = ^i%>/miA Comparison 
with Eqs. (|53"|) and ([5"4"|) shows that the transport coeffi- 
cients A and (j> predicted by the BGK model are different 
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from those obtained from Grad's method, regardless of 
the choice of the free parameter j3. 

So far, (3 has remained free. Henceforth, by following 
arguments presented in Rcfs. [45| and [46|, we will take, 
for simplicity, /3 = (l + a)/2. 

IV. SIMULATION METHODS 

As said in the Introduction, in order to assess the reli- 
ability of the previously described theoretical results and 
the existence of the LTu class, we have performed DSMC 
simulations of the Boltzmann equation and MD simu- 
lations for a granular gas of hard spheres (d — 3) [47j |. 
In the MD simulations the global solid volume fraction 
has been taken equal to 7 x 10~ 3 in order to remain in 
the dilute regime and compare with the Boltzmann re- 
sults obtained either from DSMC simulations or from 
the theoretical approaches. The gas is enclosed between 
two plates moving with velocities U± and maintained at 
temperatures T±, where the subscripts + and — denote 
upper and lower wall, respectively (see Fig. [T]). 

In our simulations we have considered N = 2 x 10 5 par- 
ticles (DSMC) and N ~ 10M0 5 particles (MD). When a 
particle collides with a wall its velocity is updated follow- 
ing the rule v — > v' + [7±x. The first contribution (v') of 
the new particle velocity is due to thermal boundary con- 
dition, while the second contribution (E/±x) is due to wall 
motion. The horizontal components of v' are randomly 
drafted from a Maxwell distribution (at a temperature 
T±) whereas the normal component v' y , due to collision 
with a wall, is sampled from a Rayleigh probability dis- 
tribution: P(\v' y \) = {m\v' y \/T ± )e- mv 'y 2 / 2T± . 

In the traditional DSMC method (2{|, which we use 
here, the system is split into cells whose characteristic 
length is much smaller than the mean free path £ (that 
is, macroscopic properties do not vary significantly along 
a cell). Here we define the (local) mean free path for 
hard spheres as £ — yjT /mv^ 1 , where the (local) effec- 
tive collision frequency v is defined in Eq. (|27p. Fur- 
thermore, the time step needs to be much smaller than 
the microscopic characteristic time (inverse of the col- 
lision frequency v). The DSMC method consists of two 
steps. One is the free streaming, where the particles move 
in straight lines without inter-particle collisions. The 
boundary conditions are applied in this step. The other 
one is the collision step, in which possible particle pairs 
are randomly selected from the same cell and collision is 
accepted with a probability 0(v.y ■ CTjj)a;y/a; max , where 
v,j = Vj — Vj is the relative velocity between particles i 
and j, ay = (r< - rj)/|r 4 - |, = (47rcr 2 n)|v ?;j • <7y|, 
and w max is an upper bound of the probability of particle 
collision per unit time. 

Given the geometry of the problem, the DSMC cells 
need not be three-dimensional since only the vertical co- 
ordinate y is recorded. This is possible because colli- 
sions are sampled independently of the particle position 
within the same layer, and only relative approach veloc- 



ities Vy • (Tij are needed in the simulation (unit vectors 
CTjj are randomly generated). In our DSMC simulations 
we have taken a time step and a layer width given by 
St = 3 x lO -3 ^ 1 and Sy = 2 x 10~ 2 £, respectively, where 
(as said in Sec. Ill CP £ = y/T- /mF 1 , v being given by 
Eq. (|27]) with T-)L and n -> n. 

In contrast to the DSMC case, a three-dimensional box 
is required in the MD simulations. We have taken h x 
hx h cubes with periodic boundary conditions along the 
directions (x and z) parallel to the walls. 

In the simulation results presented in Sec. [V] dimen- 
sionless quantities are used. We choose as units of mass, 
length, and time m, £(—h/2), and ^ _1 (— h/2), respec- 
tively, once the steady state has been reached. As said 
before, we take the condition T_ < T + . Thus, and with 
our choice of units, the reduced quantity Aj \JmT(—h/2) 
[cf. Eq. (HI])] will represent the maximum value across the 
system of the local thermal Knudscn number e [cf. Eq. 
(|25|)]. In other words, in our work the slope of the ther- 
mal Knudsen number e(y) is always positive [36jj. 

The separation between the plates has typically been 
set h si 5-20 and we have considered a wall temperature 
difference in the range AT = T+/T_-l = 0-20. Since, as 
will be seen below [cf. Fig.[9ta)], the values of the reduced 
shear rate arc smaller than 1 for a > 0.5, the above values 
of h and AT imply that the Reynolds number [cf. Eq. 
([3H| ] is always smaller than about 400. For this range 
of Re the flow is expected to remain laminar and this is 
confirmed by our simulations. 

We store instantaneous values of the relevant hydrody- 
namic quantities iterativcly at runtime, for further pro- 
cessing after the simulation run. 

With respect to the processing of the steady state hy- 
drodynamic properties, we perform two types of averages: 
one in space and the second one in time. The first one 
is performed over a number of contiguous cells, forming 
a statistical spatial bin, whose size must not be larger 
than the typical scale over which hydrodynamic fields 
vary (48[. Since this scale depends on the applied gradi- 
ents (wall temperature difference and applied shear for 
a system with a given height h), this statistical bin size 
needs to be adjusted for each simulation. We have ob- 
served that a bin adjustment of Ay w 0.1Kn _1 £ (where 
the Knudsen number is Kn = a) is enough for preserv- 
ing all properties of hydrodynamic profiles [36l [48j . The 
other averaging is performed, in each spatial bin, over 
values at different times of the same steady states. This 
double averaging is very convenient since it allows us to 
obtain very smooth hydrodynamic steady profiles, even 
if the system is not large. This is especially useful in the 
case of DSMC, where thermal fluctuations may result in 
too noisy profiles for small systems [29|, |49| . 
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FIG. 4: (Color online) Temperature vs flow velocity, T(u x ), 
as obtained from DSMC simulations at t — 45^ _1 (solid sym- 
bols, transient state) and t > 800^ _1 (open symbols, steady 
state). In these graphs AT = 4 and (a) a = 0.99 and (b) 
a = 0.4. 
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FIG. 5: (Color online) Heat flux profiles (a) q x (y) and (b) 
q y (y) as obtained from DSMC simulations at t = 45Z/ -1 (solid 
symbols, transient state) and t > 800i? _1 (open symbols, 
steady state). In these graphs AT = 4 and (a) a = 0.99 
and (b) a = 0.4. 



V. RESULTS 
A. Transient regime 

We analyze here the transition to steady LTu states 
from DSMC and MD simulation data, starting from an 
initial equilibrium distribution at T = T_. We have 
found that in general the duration of this transition to the 
steady state becomes substantially longer as inelasticity 
increases. 

Figures|4]and[5]show T{u x ) and q x ,y(y) profiles, respec- 
tively, from DSMC data for transient states at t = 45£> _1 
and steady states (t > SOOi^ 1 ) for a = 0.99 and 0.4. In 
these cases h ~ 16, as indicated by the horizontal axis of 
Fig. [5] It is apparent that, at a given common time, the 
deviations from the steady LTu profiles are weaker for 
a = 0.99 (quasi-elastic gas) than for a = 0.4 (strongly 
inelastic gas). In any case, we have seen that over the 
range of a at which we perform the simulations (a = 0.3- 
1.0), time values of about t = 250P~ 1 always yield fully 
developed steady LTu flows. This happens also for MD 
simulations, as shown in Fig. [51 where we can see results 



for temperature, heat flux, and pressure for a granular 
gas with a = 0.85. The degree of approach to the steady 
state is perhaps a little slower but, in any case, we have 
observed that the system has already reached the steady 
state at t = 250D~ 1 . 

Figures EH1 show that both DSMC and MD confirm 
the existence, in the steady state, of Couette flows with 
practically linear T{u x ) profile, uniform heat flux, and 
uniform pressure. 



B. Identification of LTu flows 

In order to identify the LTu flows, we have proceeded 
analogously to a previous work 



For each simulation 
series, we fix AT and the applied shear ([/+ — U-)/h. 
Once the steady state is reached, we monitor the para- 
metric plot of temperature versus flow field, T(u x ), look- 
ing for the typical linear profiles of the LTu steady states 
in the bulk region, that is, outside the boundary lay- 
ers. More specifically, since we observed that T(u x ) never 
shows inflection points in the bulk (in accordance with 
theory [36|]), we check the sign of the T(u x ) profile curva- 
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at both T(u x ) and q x ,y- We take as the final LTu flow 
the simulation which best approaches the conditions of 
both linear T(u x ) and constant q XtV (q z is always zero 
in our geometry). Put in other words, we find the LTu 
flows by crossing vertically (in the shearing axis direc- 
tion) the surface in Fig. [21 until getting the right value 
of the reduced shear rate (ath)- 

The degree of approach to the properties of theoretical 
LTu states that we obtained in the simulations is rather 
good. For illustration on this, we show MD simulation 
results in Fig. [7j where one can see how the transition be- 
tween states above and below the surface in Fig.[2]occurs. 
Figure Efa) shows the results for T(u x ) profiles, whereas 
Fig. [7jh) shows the corresponding results for q y (y) pro- 
files. We have found that heat flux profiles are more sen- 
sitive to a departure from the LTu surface and for this 
reason we usually proceed as described above: we first 
search for an almost linear T(u x ) profile and then we 
fine tune the LTu state by searching the flattest heat flux 
profiles for a shear rate around the first selected value. 
Compared to results from DSMC simulations (see Fig. 
2(a) in Ref. [IH ) we see that boundary layer effects on 
heat flux profiles are stronger in MD simulations. Also, 
this effect is more noticeable next to the higher tempera- 
ture wall. It is also to be noticed that the sign of d y q y (y) 
in the bulk domain changes from positive for a > a t h to 
negative for a < a t h- This agrees with the interpreta- 
tion that viscous heating carries kinetic energy toward 
the hotter wall whereas inelastic cooling tends to remove 
it from there [3(| ■ Once this effect of inelastic cooling is 
sufficiently compensated by viscous heating, we can see 
the traditional trend of heat flux profiles for clastic gases 
between two walls at different temperatures (that is, heat 
flux is directed toward the colder wall [29j]). 

In Fig. [H we show LTu heat flux profiles from DSMC 
data. It is observed that, at a given wall temperature 
difference AT, the impact of a on q y is rather weak. On 
the other hand, at a given value of a, the magnitude of 
q y is approximately proportional to AT. Although not 
shown, we have also found that the influence of a on q x 
is much stronger than in the case of q y . 



FIG. 6: (Color online) Plots of (a) T(u x ), (b) q v (y), and (c) 
p(y) as obtained from MD simulations at t = 6Qv~ l (solid 
symbols, transient state) and t > &00v~ (open symbols, sta- 
tionary state). In these graphs AT = 4 and a = 0.85. 



ture. If the sign is positive, that means that cooling still 
overcomes viscous heating. Thus, we need still increase 
the applied shear for the next simulation (while keeping 
constant AT), in the search for a T(u x ) profile with zero 
curvature. The process is repeated iteratively until we 
observe a change of sign in the curvature of T. Then, 
we look for the LTu state between the consecutive values 
for which the change of sign in the curvature is observed, 
by taking smaller changes of applied shear and looking 



C. Generalized transport coefficients 

In a recent work [HI, we introduced the method of 
measurement of the generalized transport coefficients of 
the LTu class defined by Eqs. ([SI]), JjSJ), ([34]) . and |35 ]) . 
We have confirmed by simulations that the values of these 
reduced coefficients only depend on the value of the co- 
efficient of normal restitution a. 

Figure [9] presents the simulation data for the reduced 
shear rate a, the reduced shear viscosity rj* , and the re- 
duced directional temperatures 9i as functions of the co- 
efficient of normal restitution a. The figure also includes 
the theoretical predictions obtained from Grad's method 

[cf. Eqs. jsojum and (E3 with A s ivcn b y E q- ©] 

and from the BGK-like kinetic model [cf. Eqs. Eqs. (j5t7)) - 
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FIG. 7: (Color online) Transition to LTu profiles for MD series 
with varying wall shearing at ot = 0.6. Solid symbols corre- 
spond to non-LTu states: (a) for a = 0.87a t h and (■) for 
a = 1.25a t h). Open squares (□) stand for the LTu stationary 
profile. 
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FIG. 8: (Color online) Heat flux profiles q y (y) from DSMC 
data. In panel (a) AT = 5 and a = 0.4 (□), a = 0.7 (■), 
and a = 0.99 (a). In panel (b) a = 0.7 and AT = 5 (A), 
AT = 10 (□), and AT = 15 (■). 



(5J]) and (55]) with the replacement $ x -> (1 + a)/2]. It 
can be observed a consistent agreement between DSMC 
and MD data. Moreover, the theoretical results compare 
quite well with computer simulations, the BGK results 
slightly improving the results obtained from Grad's ap- 
proximation. 

Regarding the transport coefficients characterizing the 
heat flux, Fig. [10] compares computer simulation results 
(DSMC and MD) with Grad's [cf. Eqs. (53]) and (5§ 
with & given by Eq. @U)] and BGK [cf. Eqs. (7SJ) and 
([75]) ] theoretical predictions. It is apparent that the gen- 
eralized thermal conductivity A* exhibits a weak depen- 
dence on a, in agreement with Fig. [SJa). On the other 
hand, the cross coefficient <f>* , which vanishes in the elas- 
tic limit, starts growing rapidly with increasing inelastic- 
ity, and then presents a much more moderate dependence 
on a for large inelasticities. In particular, (jf becomes 
larger than A* for a < 0.9, what represents a strong non- 
Newtonian effect. Interestingly, these features are very 
well captured by the simple Grad approximation, while 
the BGK approach only agrees at a qualitative level. The 
contrast between the good performance of the BGK pre- 
dictions for the rheological properties seen in Fig. [S] and 



the quantitative disagreement found in Fig.[TU]is in part 
due to the fact that the BGK model only possesses a 
free parameter (/?) to make contact with the Boltzmann 
equation. 



VI. CONCLUDING REMARKS 

We have presented in this work an extensive study of 
a class of granular flows recently reported @- We refer 
to this class of flows as 'LTu,' due to the linearity of 
T(u x ) profiles. Our study has been both theoretical and 
computational. In the theory part, we have presented 
results from two different approaches: Grad's moment 
method and a BGK-type kinetic model used previously in 
other granular flow problems and now applied specifically 
to the LTu flows. In the computational part, we have 
presented results also from two different methods: the 
DSMC method of the Boltzmann equation of the inelastic 
gas and MD simulations of a dilute gas. 

The objective of the paper has been two-fold. First, we 
have confirmed by computer simulations the existence of 
LTu flows in the bulk domain under strongly inelastic 
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FIG. 9: (Color online) Plot of (a) a(a) and 77* (a) and (b) 
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correspond to Grad's method and BGK model, respectively. 
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by a careful fine-tuning of the shear rate applied by the 
walls, it is possible to reach steady states with a uniform 
heat flux and a linear parametric plot of T vs u x . Sec- 
ond, we have assessed the theoretical predictions derived 
from two different approaches (Grad's moment method 
and BGK-type kinetic model) for the generalized non- 
Newtonian transport coefficients. 

The agreement for the reduced shear rate, rheological 
properties, and transport coefficients between the DSMC 
and MD simulation methods is very good, as shown in 
Figs. |9] and [TO] Also, the evolution to stationary states 
and other properties of the hydrodynamics of the LTu 
class are found to be remarkably similar for both DSMC 
and MD. Regarding the reliability of both theoretical so- 
lutions, we have observed that they are excellent for the 
rheological properties [cf. Fig. [9]. On the other hand, 
in the case of the heat flux coefficients, the quantitative 
agreement with simulation is only good for Grad's mo- 
ment method. This good performance of Grad's method 
has been also observed in the case of granular binary mix- 
tures under simple shear flow [5l|, [5S|. Nevertheless, the 
good behavior of Grad's 13-moment method does not 
extend to cases where the heat flux is not uniform, as 
happens in the Couette flow for ordinary gases [3lj |. 

As it is customary in fluid mechanics, the importance 
of describing entire classes of flows with clearly identi- 
fiable hydrodynamic properties (rather than describing 
specific properties of a given flow in a case-by-case ba- 
sis) cannot be overemphasized. In this sense, we have 
shown here that the LTu flows are characterized by a set 
of interesting properties that can be useful as a refer- 
ence point for experimental studies on granular flow at 
low density. More interestingly, we show that all flows of 
the new class share, for the same a, the same Knudscn 
number associated with transport of momentum. 

To summarize, we have described in detail the proper- 
ties of a new class of flows, finding excellent agreement 
between simulation and theory. The results show that 
this class of flows encompasses at the same time flows 
of elastic and inelastic gases, what gives solid support 
to the validity of a hydrodynamic description of granu- 
lar dynamics, at least in this case and for the type of 
geometry studied in this work. 

We expect in the future to extend these results to 
other related systems, such as mixtures, inelastic rough 
spheres, or driven systems. Also, we plan to carry out 
further studies on the hydrodynamics of this type of flows 
(instabilities, pattern formation, etc.). 



FIG. 10: (Color online) Plot of A* (a) (▲, ■) and 0(a) (A, 
□) as obtained from DSMC simulations (triangles) and from 
MD simulations (squares). The solid (A*) and dashed (<f>*) 
lines correspond to Grad's method, while the dotted-dashed 
(A*) and dotted lines (</>*) stand for the BGK model. 

conditions. At a given wall temperature difference and 
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